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ABSTRACT 

We present results of more than three decades of timing measurements of the 
first known binary pulsar, PSR B1913+16. Like most other pulsars, its rota- 
tional behavior over such long time scales is significantly affected by small-scale 
irregularities not explicitly accounted for in a deterministic model. Neverthe- 
less, the physically important astrometric, spin, and orbital parameters are well 
determined and well decoupled from the timing noise. We have determined a sig- 
nificant result for proper motion, fx a = — 1.43±0.13, /is = — 0.70±0.13 mas yr _1 . 
The pulsar exhibited a small timing glitch in May 2003, with Af/f = 3.7 x 10~ n , 
and a smaller timing peculiarity in mid-1992. A relativistic solution for orbital 
parameters yields improved mass estimates for the pulsar and its companion, 
mi = 1.4398 ± 0.0002 M and m 2 = 1.3886 ± 0.0002 M . The system's orbital 
period has been decreasing at a rate 0.997 ± 0.002 times that predicted as a re- 
sult of gravitational radiation damping in general relativity. As we have shown 
before, this result provides conclusive evidence for the existence of gravitational 
radiation as predicted by Einstein's theory. 
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Subject headings: binaries: close — gravitation — gravitational waves — pul- 
sars: individual (PSR B1913+16) — stars: fundamental parameters — stars: 
kinematics 



Introduction 



Pulsar PSR B1913+16 (PSR J1915+1606) was the first binary pulsar to be discovered 
(Hulse & Taylor 1975). Its 7.75-hr binary period and 300 km s -1 projected orbital velocity 



suggested that it would exhibit a rich set of potentially measurable relativistic effects, and 
that has turned out to be the case. Subsequent studies have shown that the pulsar's compan- 
ion must also be a neutron star. With tidal and rotational stellar distortions eliminated as 
complicating factors for any dynamical analysis, the system is an essentially clean laboratory 
for testing relativistic gravity. In this respect, the largest remaining complications depend 
on reference-frame accelerations related to the structure and dynamics of our Galaxy. 

Over the past 35 years, we and our colleagues have timed PSR 131913+16 at Arecibo 
Observatory with steadily improving equipment and analysis techniques. Among the best 
known results are measurement of the general relativistic advance of periastron at a rate 



some 35,000 times that of Mercury in the solar system (Taylor et al. 1976); the effect of 



gravitational radiation damping, causing a measurable rate of orbital decay (Taylor et al. 



1979); and detection of changes in the pulse shape, resulting from geodetic spin precession 



(Weisberg et al. 1989). Our results continue to be fully consistent with general relativity 



and have placed strong constraints on alternative, previously viable, relativistic theories of 
gravity ( |Weisberg & Taylor|1981| |Taylor fc Weisberg|l982| p§9| |Taylor et al.[l992| |Weisberg 
fc Taylor|[2002| |Clifton fc Weisberg|[2008| ) . The purpose of this paper is to provide our latest 
analysis of the pulsar and its orbit, based on timing observations from 1974 through 2006. 



2. Observations 

Our highest quality dataset consists of 7650 five-minute integrations of the pulse pro- 
file at wavelengths near 20 cm, each subsequently reduced to a topocentric time-of- arrival 
(TOA) measurement at Arecibo Observatory. These data span the years 1981 through 2006; 
earlier measurements are consistent with these but (owing to their much larger uncertainties) 
have little effect on the final results. In 2004 we made a special effort to observe regularly 
throughout the year, so as to improve quality of the astrometric results. Parameters of the 
observing equipment and associated measurement uncertainties are presented in Table [I] 
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Further details of the Princeton Mark I observing system are given by |Taylor fc Weisberg" 



(1982), while the Mark II and III are described in Taylor & Weisberg (1989). The latest 
system, the Wideband Arecibo Pulsar Processor (WAPP), consists of four spectrometers 
optimized for pulsar work, each having 512 frequency channels across a 100 MHz bandwidth 
and synchronously accumulating the pulsar signal into 1024 phase bins across the full pulsar 
period. 



3. Analysis of TOAs 

Data were analyzed with the program TEMPO, using the JPL DE405 solar system 



ephemeris and the relativistic timing model of Damour & Deruelle (1986). The Damour- 
Deruelle approach is particularly desirable because it leaves open the question of the correct 
relativistic theory of gravity: one solves for phenomenological parameters whose precise 
relation to a gravitational theory can be investigated later. Further discussion of the Damour- 



Deruelle model and its application to binary pulsars in TEMPO can be found in Taylor Sz 



Weisberg (1989). Important quantities determined in this way include those related to the 
pulsar's celestial position, spin, and orbit. The first category includes right ascension a, 
declination 5, proper motions fi a and epoch to, pulse repetition frequency /, spindown 
rate /, and a glitch epoch and frequency discontinuity Af. Fitted orbital parameters include 
five Keplerian quantities: the projected semimajor axis of the pulsar orbit x = aj.sinz, 
eccentricity e, epoch of periastron passage T , period Pb, and longitude of periastron coo; and 
the following relativistic or "post-Keplerian" parameters: average rate of periastron advance 
(u), variations in gravitational redshift and time dilation 7, and orbital period derivative JV 
As with most other pulsars that have been timed carefully over several decades, a number of 
"nuisance parameters" must also be measured to account for unmodeled long-term timing 
irregularities. These extra fitted terms have no clear physical interpretation beyond being 
somehow related to the poorly understood structure and dynamics of the spinning neutron 
star. Their mathematical form is somewhat arbitrary; as described further below, we have 
experimented with a number of higher-order time derivatives of the pulsar rotation frequency, 
and one additional frequency discontinuity. 



3.1. Astrometric and Spin Parameters 

The astrometric and basic spin parameters determined from the full data set are listed 
in Table [2] A significant result for proper motion has been obtained for the first time. While 
the pulsar position is determined from the annual variation of TOAs as the Earth moves 
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about its orbit, proper motion is determined from a secular variation of this annual signal. 
Timing noise or other systematic effects can contaminate such measurements. We now believe 



that the previously reported value of PSR B1913+16's proper motion (Taylor & Weisberg 



1989) was biased by timing noise, a problem which was exacerbated by the concentrated but 



~biennially spaced observing campaigns frequently employed for this source. To circumvent 
these problems, we observed PSR B1913+16 several times over the course of calendar year 
2004 to achieve thorough data coverage around the Earth's orbit. By merging the 2004 data 
with observations made in 1985-1988, which also had good coverage throughout those years, 
we have now obtained a robust measurement of the pulsar's proper motion. See Section [5] 
for further analyses of the proper motion result. 

As in other pulsars, the measured spindown rate / is assumed to be the result of an 
electromagnetic braking torque on the spinning, strongly magnetized neutron star. Its value 
is deterministic and (at the level of accuracy quoted in Table 2) independent of the span 
of data over which it is fitted. In addition, the pulsar experienced a well-defined "classical" 
glitch in May 2003. The magnitude A/// = 3.7 x 10 -11 of the event is smaller than almost 
all seen in the population of normal (not recycled) pulsars, with only some of the glitches 



in PSR B0355+54 having a comparably small magnitude (Melatos et al. 2008 Chukwude 



& Urama 2010). Ours is only the second glitch to be detected in a recycled pulsar, with 



the other (the smallest known among all pulsars) occurring in the millisecond pulsar PSR 



B1821-24 in globular cluster M28 (Cognard & Backer 2004) 



The remaining timing parameters — higher-order derivatives /, /,..., and (in some fits) 
another small frequency discontinuity in mid- 1992 — were introduced to the timing solution 
in an ad hoc manner, in order to "whiten" the remaining post-fit residuals. Although we 
offer no clear or unique physical interpretation for these parameters, their combined effects 
are almost certainly a consequence of stochastic timing-noise processes in the neutron star 
interior ( Cordes fc Downs|1985 Arzoumanian et al.|1994 Urama et al.|2006 ). Unlike the case 
of the well-fitted May 2003 glitch, the values of these fitted parameters are not independent of 
the data span analyzed and cannot be expected to extrapolate the timing behavior accurately 
to future epochs. Identifying the mid-1992 behavior as a discrete event is highly uncertain, 
in part because of coarse sampling around that time. Our data can be fit almost as well 
by introducing several additional frequency derivatives instead of a second discrete event. 
Similar results were obtained by fitting multiple harmonically related sinusoids to the timing 



noise with the routine "FITWAVES" (Hobbs et al. 2004) of program TEMP02 (Hobbs et 



aL||2006 Edwards et al.||2006 ), instead of the multiple frequency derivatives described above. 
Since we are unable to converge on a unique glitch parameter solution for the mid-1992 
behavior, we do not include this event in Table [2} 
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3.2. Orbital Parameters 



Fitted orbital parameters for PSR B1913+16 are listed in Table |3j As we have empha- 



sized before (Taylor & Weisberg 1989 ), values for each of the Damour-Deruelle post-Keplerian 



parameters expected in general relativity can be expressed in terms of the Keplerian param- 
eters and the initially unknown masses of the pulsar and its companion, mi and m 2 . The 
appropriate expressions for {Co) and 7 are 



3 G 2 ' 3 c~ 2 (P fe /2vr) 



-5/3 



2.113323(2) 



(mi + m 2 ) 



1-e 2 )- 1 

2/3 

deg yr _1 



m-i + m 2 



(2/3 



7 = G 2/3 c~ 2 e (P b /2n) 1/3 m 2 (m x + 2m 2 ) (m x + m 2 )" 4/3 



0.002936679(2) 



m 2 (mi + 2m 2 ) (mi + m 2 ) 



-4/3 



AC 



2/3 



s. 



(2) 



In the second line of each equation we have substituted values for P b and e from Table [3j 
and used the constants GM Q /c 3 = 4.925490947 x 10~ 6 s and 1 Julian yr = 86400 x 365.25 s. 
The figures in parentheses represent uncertainties in the last quoted digit, determined by 
propagating the uncertainties listed in Table [3} In each case, the uncertainties are dominated 
by the experimental uncertainty in orbital eccentricity, e. 

Eq. (1) may be solved for the total mass of the PSR B1913+16 system, yielding M = 
m i + m 2 — 2.828378 ± 0.000007 M . The additional constraint privided by Eq. (2) permits 
a solution for each star's mass individually, mi = 1.4398 ± 0.0002 M and m 2 = 1.3886 ± 
0.0002 M . As far as we know, these are the most accurately determined stellar masses 
outside the solar system. It is interesting to note that since the value of Newton's constant 
G is known to a fractional accuracy of only 1 x 10~ 4 , M can be expressed more accurately 
in solar masses than in grams. 



3.3. Gravitational Radiation Damping 



According to general relativity a binary star system should radiate energy in the form 
of gravitational waves. Peters and Matthews (1963) showed that the resulting rate of change 
in orbital period should be 

192 ttg 5 / 3 fp b y 5/3 

5^5 \2k) 

_ 12 mi m 2 (mi + m 2 ) -1 / 3 



pGR 



1 



73 



37 



— e* H e 

24 96 



1 



= 2x-7/2 



mi m 2 (mi + m 2 ) 1 ^ 3 



-1.699451 8 x 10" 



M, 



5/3 







(3) 
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Inserting values obtained for m x and m 2 and propagating uncertainties appropriately, we 
obtain the general relativistic predicted value 



pGR 



-2.402531 ± 0.000014 x 10 



-12 



(4) 



Equations (3) and (4) apply in the orbiting system's reference frame. Relative acceler- 
ation of that frame with respect to the solar system barycenter will cause a small additional 
contribution to the observed p,. Damour & Taylor (1991) presented a detailed discussion 



of this effect and other possible contributions to P&. Recent progress in determining the 
galactic-structure parameters allows us to update the relevant quantities and compute a new 
value for the kinematic correction to P&. Using R Q = 8.4 ± 0.6 kpc for the distance to the 
galactic center and Go = 254 ± 16 km s" 1 for the circular velocity of the local standard of 
rest (IGhez et al.||2008| IGillessen et al. 1120091 iReid et al. 1120091), and d = 9.9 ± 3.1 kpc for the 



pulsar distance (Weisberg et al. 2008), we obtain the kinematic contribution, APb iga i: 



AP b 



,gal 



-0.027 ±0.005 x 10 



-12 



Thus, we find the ratio of observed to predicted rate of orbital period decay to be 



AP b . 



gal 



pGR 



0.997 ±0.002. 



(5) 



(6) 



Agreement between the observed orbital decay and the general relativistic prediction is 
illustrated in Fig. |2j which shows how excess orbital phase (relative to an unchanging orbit) 
has accumulated since the pulsar's discovery in 1974. We note that the overall experimental 
uncertainty embodied in Eq. (6) is now dominated by uncertainties in the galactic parameters 
and pulsar distance, not the pulsar timing measurements. Even better agreement between 
observed and expected values of P& would be obtained if the true value of R or d were 
slightly smaller, or Bo slightly larger. For example, observed and expected values agree if 
d = 6.9 kpc, which is within the Weisberg et al. (2008) error envelope. It will be interesting 



to see whether improved future estimates of these quantities will show one or more of these 
conditions to be true. 



Other Relativistic Effects 



Two other relativistic phenomena are potentially measurable in the PSR B1913±16 
system: geodetic precession and gravitational propagation delay. Spin-orbit coupling should 



Barker & O'Connell 


1975a 


shape changes. 


Weisberg et 



al. ( 1989 ) first detected such changes, which were observed and modeled further by Kramer 



(1998). Weisberg & Taylor (2002) and Clifton & Weisberg (2008) found that the pulsar 



beam is elongated in the latitude direction and becomes wider in longitude with increasing 
distance from the beam axis in latitude. These models suggest that in the next decade or so, 
precession may move the beam far enough that the pulsar will become unobservable from 
Earth for some decades, before eventually returning to view. 



The excess propagation delay (Shapiro 1964 ) caused by passage of pulsar signals through 



the curved spacetime of the companion is largest at the pulsar's superior conjunction. The 
maximum amplitude varies with time because the impact parameter at superior conjunction 
depends strongly on the current value of u. In this respect the orbital geometry was partic- 
ularly unfavorable in the mid-1990s (see Damour & Taylor 1992), but in coming years the 



propagation delay should start to become observable. Damour & Deruelle (1986) charac- 
terize the measurable quantities as range r = (Grri2/c 3 ) and shape s = sini, of the Shapiro 
delay, where i is the orbital inclination. As orbital precession carries our line of sight deeper 
into the companion's gravitational well, future observations should permit the robust mea- 
surement of these two parameters, and hence two additional tests of relativistic theories of 



gravity (Damour 2009 Esposito-Farese 2009). 



5. Systemic Velocity 



Our pulsar proper motion measurement (Section 3.1), combined with the distance es 



timate discussed in Section 3.3, corresponds to a transverse velocity (with respect to the 
solar system barycenter) of 75 kms -1 with a galactic position angle of 306°; i.e., directed 36 
degrees above the galactic plane. The ~ 30% distance uncertainty places similar limits on 
velocity accuracies. 

We can now estimate two components of the pulsar systemic velocity in its own standard 
of rest by combining the measured pulsar transverse velocity and distance, the solar motion 



with respect to our Local Standard of Rest (Schonrich et al. 2010) , and galactic quanti- 



ties Ro and Go- The third component of motion, which is inaccessible via proper motion 
measurements, lies close to the direction of Galactic rotation at the pulsar's position. 

The pulsar's galactic planar and polar velocity components relative to its standard of 
rest are 247 km/s almost directly away from the galactic center and 51 km/s toward the 
galactic North Pole, respectively. (This is significantly larger than the measured velocity in 
the solar system barycenter frame because the pulsar's standard of rest velocity fortuitously 
cancels much of the pulsar's peculiar velocity with respect to it.) The systemic velocity 
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of B 19 13+16 is significantly larger than other well-measured double neutron star binary 
system velocities, including the J0737-3037 (transverse velocity 10 km/s; Stairs et al. 2006), 
J1518+4905 (transverse velocity 25 km/s; Janssen et al. 2008), and B1534+12 (transverse 
velocity 122 km/s; Thorsett et al. 2005) systems. 

6. Conclusions 

We have analyzed the full set of Arecibo timing data on pulsar B1913+16 to derive 
the best values of all measurable quantities. A significant proper motion has finally been 
determined. A small glitch was observed in the pulsar's timing behavior, the second known 
glitch in the population of recycled pulsars. The measured rate of orbital period decay 
continues to be almost precisely the value predicted by general relativity, providing conclusive 
evidence for the existence of gravitational radiation. Uncertainties in galactic accelerations 
now dominate the error budget in P b , and are likely to do so until the pulsar distance can 
be measured more accurately. We expect that the Shapiro gravitational propagation delay 
will yield additional tests of relativistic gravity within a few more years. 

The three authors gratefully acknowledge financial support from the US National Sci- 
ence Foundation. Arecibo Observatory is operated by Cornell University under cooperative 
agreement with the NSF. We thank Joseph Swiggum for assistance with analyses of glitches 
in the pulsar population; and C . M. Ewers^J A. de la Fuente, J.T. Green, and Z. Pei for 
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Table 1. Parameters of Observing Systems 



Name 


Dates 


Total 


Frequency 


Time 


TOA 


Number < 






Bandwidth 


Channels 


Resolution 


Uncertainty 


TOAs 






(MHz) 




(H 


(/is) 




Mark I 


1981 - 1984 


16 


64 


125 


20 


1719 


Mark II . . . 


1985 - 1989 


8 


32 


125 


31 


1015 


Mark III. . . 


1988 - 2003 


40 


32 


640 


16 


2473 


WAPP a . . . 


2003 - 2006 


100 


512 


58 


14 


2443 



a Three or four WAPPS were simultaneously employed in nonoverlapping frequency bands. 
The listed numbers refer to a single WAPP. 



Table 2. Astrometric and Spin Parameters 



Parameter 


Value a 


t (MJD) b 


52984.0 


a (J2000) 


19H5 m 27;99928(9) 


5 (J2000) 


16°06'2773871(13) 


H a (mas yr" 1 ). . . 


-1.43(13) 


(x s (mas yr" 1 ) 


-0.70(13) 


/(s- 1 ) 


16.94053778563(15) 


/(s- 2 ) 


-2.4761(9) xl0~ 15 


Glitch epoch (MJD) 


52770(20) 


A/ (s- 1 ) 


6.2(2) xl0~ 10 



a Figures in parentheses represent estimated uncertainties in the last quoted digit. The 
estimated uncertainties range from (3 — 10) x the formal fitted uncertainties, in order to also 
reflect variations resulting from different assumptions regarding timing noise, etc. 

b This quantity is the epoch of the next six measurements tabulated here. 
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Table 3. Orbital Parameters 



Parameter 


Value a 


T (MJD) 


52144.90097841(4) 


x = ai sin 2 (s). . . 


2.341782(3) 


e 


0.6171334(5) 


n (<i! 


0.322997448911(4) 


w (deg) 


292.54472(6) 


{Co) (deg / yr) . . . 


. 4.226598(5) 


7 (ms) 


4.2992(8) 


n 


-2.423(1) xlO- 12 



a Figures in parentheses represent estimated uncertainties in the last quoted digit. The 
estimated uncertainties range from (2 — 6) x the formal fitted uncertainties, in order to also 
reflect variations resulting from different assumptions regarding timing noise, etc. 
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Fig. 1. — Timing residuals for PSR B1913+16. (a) Residuals from a fit for data before 
mid-1992. The glitch in May 2003 can be recognized by a distinct change in slope of the 
residuals versus time. The apparent change in mid-1992 is much smaller and may or may 
not involve a discrete event, (b) Residuals from a fit of all data, holding astrometric and 
orbital parameters fixed at the values in Tables [2] and [3] fitting for pulsar frequency and 
spin-down rate, / and /; and not allowing for higher order frequency derivatives or glitches. 
The glitch in May 2003 is evident as a sharp discontinuity, (c) Residuals from the full timing 
fit, including higher order frequency derivatives and the glitch. 
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1975 1980 1985 1990 1995 2000 2005 

Year 



Fig. 2. — Orbital decay caused by the loss of energy by gravitational radiation. The parabola 
depicts the expected shift of periastron time relative to an unchanging orbit, according to 
general relativity. Data points represent our measurements, with error bars mostly too small 
to see. 



